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In this paper we apply to gravitational waves (GW) from the merger phase of binary systems a 
recently derived frequentist methodology to calculate analytically the error for a maximum likelihood 
estimate (MLE) of physical parameters. We use expansions of the covariance and the bias of a 
MLE estimate in terms of inverse powers of SNRs where the square root of the first order in the 
covariance expansion is the Cramer Rao Lower Bound (CRLB). We evaluate the expansions, for 
the first time, for GW signals in noises of GW interferometers. The examples are limited to a 
single, optimally oriented, interferometer. We also compare the error estimates using the first two 
orders of the expansions with existing numerical Monte Carlo simulations. The first two orders of 
the covariance allows to get error predictions closer to what is observed in numerical simulations 
than the CRLB. The methodology also predicts a necessary SNR to approximate the error with the 
CRLB and provides new insight on the relationship between waveform properties, SNR, dimension 
of the parameter space and estimation errors. For example the timing match filtering can achieve 
the CRLB only if the SNR is larger than the Kurtosis of the gravitational wave spectrum and the 
necessary SNR is much larger if other physical parameters are also unknown. 



I. INTRODUCTION 

The ground-based gravitational waves detectors 
LIGO, Virgo, and GEO 600 (P, [2], [3]) are rapidly 
improving in sensitivity. By 2015, advanced versions 
of these detectors should be taking data with a design 
sensitivity approximately 10 times greater than the 
previous generation, and the probed volume will grow by 
a factor of about a thousand. Such improvements in de- 
tector sensitivity mean that the first gravitational-wave 
signature of a compact-binary coalescence (CBC) could 
be detected in the next few years (see for example [4]). 
Among the expected signals, a special role is covered 
by inspiralling compact binaries. This follows from the 
ability to model the phase and amplitude of the signals 
quite accurately and consequently to maximize the 
signal-to-noise ratio (SNR) by using matched filtering 
techniques. Matched filters also provide a maxim 
likelihood estimation of the waveform parameters such 
as component masses or time of coalescence. The choice 
of the MLEs as reference estimators is also motivated 
by the fact that if an unbiased estimator that attain the 
CRLB exists, is the MLE154 . 

The existing GW frequentist literature ([I5]..[55]l 
evaluates the MLE accuracy in two ways: (a) analyt- 
ically by calculating the so called Fisher information 
matrix (FIM) or equivalently the Cramer Rao Lower 
Bound (CRLB) which is the square root of the diagonal 
elements of the inverse FIM, and (b) numerically by 
performing Monte Carlo simulations. 
The FIM was derived analytically in ( [33] , [SD] , [U ) 
using Newtonian waveforms, extended to second-order 
post-Newtonian ( [bQ_ , [_55 ) and recently revisited up to 
3.5PN (@6],[47]). In (^),([47]) the authors calculate 
the CRLB for the three standard binary systems (NNS, 



NBH, BBH), and show how the errors change when the 
successive different PN orders are taken into account. 
They consider initial LIGO, advanced LIGO and Virgo 
noises. They also considers PN corrections to the 
amplitude. 

Monte Carlo simulations were performed, for example 
in ([IH], [H]), for the lowest PN orders, where it is also 
suggested that the inclusion of the higher PN orders 
would be computationally expensive. More recent Monte 
Carlo simulations with 3.5 PN waveforms are described 
in (^ ). 

We did not try to compare the uncertainties derived 
here to other existing papers (especially from the '90) 
since different parameter sets, noise spectra (even for the 
same antenna) and PN terms were used. For example 
in ([S^) a comparison between the CRLB and other 
bounds is done for a waveform at the 0th PN order. This 
work also uses different conventions on the waveform 
spectrum than more recent literature. In ([S3]) phasing 
is extended to the 1.5PN order. The spin parameters 
are taken in account but the noise spectrum for LIGO 
is different than the currently used design noise curves. 
In ([Sg), ([5S]), (gg) the 2PN wave was used. In the 
work ([55') interesting observations are made about 
the fluctuation of the parameters variance with the PN 
orders, analyzing both the case of spin and spinless 
systems. The fluctuations of the variance in the spinless 
case is also stressed in ([^). 

The CRLB is a convenient tool to approximate the 
accuracies in large SNRs and to obtain error bounds 
for unbiased estimators. Unfortunately, for low SNRs 
(below 20) where the first detections might emerge, 
the CRLB can grossly underestimate the errors 
[S2] , [33] , [3S] , [53] , [3H] • The reason is that with non linear 
data models and (or) with non Gaussian noises the 
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CRLB depends only on the curvature of the hkelihood 
function around the true value of the parameter. 

In this paper we apply a recently derived analyti- 
cal tool to better predict an MLE accuracy and to 
establish necessary conditions on the signal-to-noise 
ratio (SNR) for the MLE error to attain the CRLB. Ex- 
plicitly, within the frequentist framework, for arbitrary 
probability distribution of the noise, expansions of the 
bias and the covariance of a MLE in inverse powers of the 
SNR are discussed. The first order of the expansion of 
the variance is the inverse of the FIM. By requiring that 
the second order covariance is smaller, or much smaller, 
than the first order this approach predicts a necessary 
SNR to approximate the error with the CRLB. The 
higher order of the expansions are determinant in the 
low SNR regime where the inverse FIM underestimates 
the error. We compared our the errors computed using 
the first two orders of the expansions to the Monte Carlo 
simulations in ([5l]). We observed the first two orders 
of the expansions get error predictions closer than the 
CRLB to what is observed in the numerical simulations. 
In ([Sni) the simulations are related to the FIM to 
estabish ranges of SNR where the CRLB describe the 
error. Our expansions predict the same SNR range of 
validity for the CRLB. 

The expansions are sensitive to the side lobes of the 
likelihood function because they make use of higher order 
derivatives than the second one (which is only sensitive 
to the curvature of the main lobe). The methodology 
also provides new insight on the relationship between 
waveform properties, SNR, dimension of the parameter 
space and estimation errors. For example the timing 
match filtering accuracy achieves the CRLB only if the 
SNR is larger than the Kurtosis of the gravitational 
wave spectrum and the necessary SNR is much larger 
if the other physical parameters are unknown. More 
specifically the MLE of the arrival time for NS-NS 
binary signals might require an SNR equal to 2 with 
the time as the only parameter or 15 when all the other 
parameters are unknown. These results are important 
to understand the domain of validity of recent papers 
like fS7| that defines 90% confidence regions in direction 
reconstruction with time triangulation. The regions 
discussed in [ST] for SNR smaller than 10 are based on 
timing MLEs, with the arrival time being the only un- 
known parameter, and the time uncertainties quantified 
by the CRLB. 

We also note that |34i, using a formalism introduced in 
[5S],[33], describes a consistency criterion, different from 
the condition derived in this paper, for the validity of the 
CRLB that, if applied to a 4pp compact binary signal 
computed with a 2PN expansion and mi — m2 — IOMq, 
requires an SNR of at least 10. At the time of the 
writing of this paper we established with M.Vallisneri 
that the equation ( |3.19[ ) of this paper becomes, in 
the one parameter case and colored Gaussian noise, 
equivalent to equation (60) in ([M]) or (A35) in ([3S]). A 
comparison for the Gaussian noise and multi parameter 



case is object of current work, while a comparison 
for arbitrary noise is not possible because ([34j) and 
([3S]) use from the beginning of their derivations Gaus- 
sian noises. The explicit calculations shown here for 
different GWs are also not performed in ((34j) and f(35j). 

In section |TT] we present the explicit expressions of 
the expansions of the bias and the covariance matrix 
for arbitrary noise and size of the parameter space. 
In section |III| we explain how the expansion can be 
evaluated for signals in additive colored Gaussian noise. 
In section |IV| we describe the post-Newtonian inspiral 
waveform used for the examples, describe the parameter 
space and the initial and advanced LIGO noises. In 
section [V| we study the one-dimensional parameter space 
results when only one parameter at a time is considered 
unknown. In section VI We present the results for full 



parameter space with the initial and advanced LIGO 
noises. We also compare our results with published 
results from Monte Carlo simulations. In section IVIII 
we present some conclusions and in the appendix we 
describe the derivation of the expansions as well as the 
relationship of this method with the statistics literature. 



II. EXPANSIONS FOR THE BIAS AND 
COVARIANCE MATRIX OF A FREQUENTIST 
MLE IN ARBITRARY NOISE 

In this section we present the first two orders of the 
expansions in inverse powers of the SNR for the bias 
and the covariance matrix. The details of the derivation 
are provided in appendix A. Given an experimental data 
vector X = {xi, ..,xn}, where N is the dimension, we 
assume that the data are described by a probability den- 
sity P{x, '&) that depends on a D-dimensional parameter 
vector d = {t^i, .., i?/?}. According to ([E]), we suppose 
that the MLE = .., i^u} of is given by a station- 
ary point of the likelihood function l{x,'&) = ln{P{x,^)) 
with respect to the components of t? 

/.fel)=^^l,.^ = r = l,..,7^. (2.1) 



If we introduce the notations 
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where —Vab is the Fisher information matrix iab — 
—Vab = —E[lab] = E[lalb] {E[.] IS the expected value), 
the first two orders of the bias for the MLE of the r 
component of the parameter vector § are given by 
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were we assumed the Einstein convention to sum over 
repeated indices. For the covariance matrix the first or- 



der is the inverse of the Fisher information matrix while 
the second order is given in by (for simpUcity we provide 
the diagonal terms): 
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III. EXPANSIONS FOR SIGNALS IN 
ADDITIVE COLORED GAUSSIAN NOISE 

For this analysis we assume that the output of an op- 
timally oriented GW interferometer has the form: 



where h(t, 9) is the signal, which depends on the param- 
eters vector 0, and w{t) a stationary Gaussian noise with 
zero mean. The probability distribution can be written: 



x{t) = h{t,e)+w{t) 



(3.1) 



p(a;)aexp|-^ J [x{t) - h{t,e)]Cl{t - ti)[x{ti) - h{ti,e)]dtdti^ 



(3.2) 



The first and second derivative of the log-likelihood give: 
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^ab = E[lalb] = -E[u] = J Kit,e)n{t-ti)hb{h,e)dtdh 

dfdf'ha{f)ht{f)n{-f,-f) 



(3.5) 
(3.6) 



It's easy to verify that 
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where we have introduced, S'/i(/),the one sided power 
spectral density defined as the Fourier transform of the 
noise auto-correlation: 



R{t) = E[n{t + T)n{T)]. 
5„(/) = J dte-^^'f'R{t). 



(3.8) 
(3.9) 



Notice that the sign convention in (3.7) follows from the 



implicit assumption that R(t) is the Fourier transform of 
E[n{f)n{f')]. In the literature another convention with 
the minus sign is also found corresponding to defining 
R(t) as the Fourier transform of E[n{f)n* {f')]. Using 
the relation h{—f ) = h{f)*, we can finally write the FIM: 



lab^E[lalb]^{haif),h{f)) 



(3.10) 



where ha{f) are the derivatives of the Fourier transform 
of the signal with respect to the a-th parameter. We have 
introduced a mean in the frequency space: 



{u{f) , v{f)) ^ m 



u{f)v{f)* 
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where the range of integration depends on the antenna 
properties and the theoretical model for the binary sys- 
tem. 

The SNR corresponds to the optimal filter: 

We can express in the same way all the quantities we 
need in order to calculate the second order variance, like 
scalar products of h{f) derivatives. 

Va.b = -Vab = iab = {K , h) (3-12) 
'>Jab,c^ {hab , he) (3.13) 
Vabc = ~-{hab , he) ~ {hac , h) - {he , ha) (3.14) 

Vab . cd = {hab , hcd) + {ha , h) {h^ , hd) (3.15) 
Vabc,d = {habc , hd) (3.16) 

Vabcd = — {hab , hcd) — {hac , hd) — {had , he) — 
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- {hacd , h) - {hcd , ha) (3.17) 

Vab,c,d^ {ha, h) {h , hd) = iab led (3.18) 

If one uses these functions, the form for the second 



(3.11) order variance, eq (2.4) can be further simplified: 
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(3.19) 



IV. INSPIRAL PHASE WAVEFORM FOR 
BINARY SYSTEMS 

We apply the general results of the theory to the case 
of a binary system in the inspiral phase. Starting from 
the the 3.5PN phasing formula ([43]) we write the Fourier 
transform of the chirp signal: 



h{t) = a{t) 



(4.1) 



where </>(/) is the implicit solution of the 3.5PN phasing 
formula, using the stationary phase approximation (SPA) 
(m, iSI)- The final result is: 



hif) - Af- 
The phase is given by: 
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(4.2) 
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where t and are the arrival time and the arrival phase. 
The function v can be defined either in terms of the total 
mass of the binary system, M = mi + m2, or in terms of 
the chirp mass, A4 = AI: 
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where 77 is the symmetrized mass ratio 
mim2 



The amphtude is a function of the chirp mass, the effec- 
tive distance, and the orientation of the source: 

A oc M'^Q{angles)/D. 

The coefficients a^'s with k = 0..N (the meaning of each 
terms being the | PN contribution) are given by: 

ai = 0, 
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where A ~ —0.6451, 6 ~ —1.28. 7 is the Euler constant, 
and Vlso — {t^M fiso)^ , with fiso the last stable orbit fre- 
quency for a test mass in a Scharzschild space-time of 
mass M: 



fiso = {6-2 ttM)- 



(4.4) 



which wiU be also used as upper cutoff for the integrals 

(IstTI) .. pl8|. 



Given the waveform (4.2 1 one can easily calculate the 
Fisher information matrix, and its inverse, the CRLB. 



(4.2) contains five unknown parameters, {A,t,(j),Ai,r]) 
(the total mass M could be used instead of the chirp 
mass), so that one should calculate a five dimensional 
square matrix. 

It was already observed by (eg [S3] ) that the errors 
in the distance, and consequently the amplitude A, 
are uncorrelated with the errors in other parameters, 
ie that the Fisher information is block diagonal. We 



observed that this is also the case for the waveform (4.2) 



we use here for both the FIM and the second order 
covariance matrix. We can therefore work in a four 
dimensional parameter space (t, (j), M,r]). However it is 



worth observing that this is not an obvious point since 
in general the amplitude estimation can be coupled with 
other other parameters if they enter in the signal in 
specific ways (see ch.3 of |54) ) . 

It is also worthwhile noticing that the SNR and the 
amplitude A are related like follows: 



p'^{h{f),h{f))=4A' 
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We perform the calculations using analytical forms of 
the design initial and advanced LIGO noise spectrum 
( 4.6 ) .The initial one sided power spectral density of LIGO 



can be written for / > (Shif) = 00, / < flow)- 



S,,{f)) = So (4.49a;)"^^ + 0.16x'*-^^ + 0.52 + 0.32a; 

(46) 

Where the lower frequency cutoff value is flow = 40Hz, 
x= j^, fo = 500Hz, and Sq = 9 x IQ-'^'^Hz-^ 

The Advanced LIGO one sided psd has the following 
expression, for / > flow {Sh{f) =00, / < flow)- 
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Where the lower frequency cutoff value is flow = 20Hz, 



X = jr, So — 215IIz, and 5*0 
plotted in fig. [T] 
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Figure 1. The Initial (dashed line) and Advanced (solid line) 
LIGO noise spectrum. 

We now calculate the second order covariance for the 
full four dimensional parameter space (t, 0, A^, ry). 

In order to make easier the comparison of our results 
with the literature we study a binary neutron star sys- 
tem (NS-NS), a neutron star - black hole system (NS- 
BH), and a black hole system (BH-BH). A neutron star 
is assumed to have a mass of 1.4Af0 and a black hole of 
IOA/0 
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We performed our calculations following these steps: 
(I) We give to the total mass M and the mass ratio 77 a 
numerical value depending on the binary system we con- 
sider. This makes the upper cutoff ( 4.4 ) to have a numer- 



ical value. (II) We compute analytically the derivatives of 
the wave function h[f) and use it to compute the Fisher 
information ( 3.10p and his inverse, the CRLB. (Ill) We 
calculate the w's functions (3.12) .. (3.18) and use them 



to compute the second order covariance (3.19) (IV) We 
plot for the four parameters the CRLB, the second or- 
der variance and the total error (without bias, see end of 
section W\ 



V. ONE DIMENSIONAL PARAMETER SPACE 

In this section we describe the results for the instance 
where only one of the parameters (we call it 0) is con- 
sidered unknown, while the others have known values. It 
can be shown that in this case the second order variance 



(3.19) can be written as: 



C2{e) = {I'y {Si'^ihee , hef ~ {heee , he)) (5.1) 

Let us consider the case were 9 is the arrival time t in 
the waveform (4.2). The ratio between the second order 



and the first order variance turns to be: 
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|/i(/)P is the a-th moment of 



the signal's frequency spectrum. The second order is neg- 
ligible with respect to the first if 



(5.3) 



that is, if the SNR is much larger than the Kurtosis of 
the signal's spectrum. This means that for two signals 
with the same bandwidth, the one with a more peaked 
spectrum becomes requires higher SNR to attain the 
CRLB. It must be noted that the Ka are functions of 
the total mass via the upper limit of the integral (see 
(Pl» 



It can be shown that for a BH-BH (NS-NS) system the 
second order becomes equal to the first for p = 1.32 (p = 
1.97). These values of SNR are smaller than those we 
will derive with the full parameters space. This indicates 
that it is much harder to estimate all the parameters 
at the same time. Also notice that if someone requires 
the second order to be much smaller than the first, for 
example 10% the conditions become more stringent. 



A similar calculation can be done when one of the 
other parameters is considered as unknown. For M and 
77 the same analysis does not give a result that is equally 
compact and we only show a plot of the ratio between 
the second order and the first order variance, for a fixed 
value of SNR, p = 10; These values, are presented in Fig. 
[2] for different values of rj and M and IPN waveforms. 
When rj is the unknown parameter (upper plot), the ra- 
tio becomes smaller when the total mass of the system 
increases. This is in agreement with the Monte Carlo 
simulations (see e.g. [55], and section 8.3) performed in 
the full four parameter space. The necessary SNR is not 
very sensitive to the value of rj. Similar considerations 
can be drawn when the total mass M is the unknown 
parameter (bottom plot), except that now the necessary 
SNR seems to be slightly more sensitive to the value of 
77. In both panels[2]the second order is always about 30% 
of the first order. If one works with a one dimensional 
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Figure 2. The ratio (72(77)/^'"' (top) and C2 (M)/i^'^*^ (bot- 
tom) 



parameter space more compact expressions can also be 



given for the first and second order bias, eqs. (2.2) and 
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We observed that the first two orders of the bias are, for 
all the the examples of this paper few order of magnitude 
smaller than the variance contributions. Therefore we do 
not include them in the presented results. Ongoing re- 
search on waveforms including the merger and ringdown 
phases show that the bias can become important too for 
systems with masses beyond the range considered here. 



VI. FULL PARAMETER SPACE 

We present results for the full parameter space begin- 
ning with the examples obtained using the initial LIGO 
noise (4.6) and then show the same results for the ad- 
vanced LIGO noise (4.7). In each plot we show three 
curves as a function of the SNR. The dotted one is the 



CRLB (square root of the diagonal elements of the in- 
verse of the FIM). The dashed one is the square root of 
the corresponding diagonal elements in the second order 
covariance matrix, and the continuous one is the square 
root of the sum of the diagonal elements of the of the FIM 
and of the second order covariance matrix. For all the 
cases analyzed in this paper the bias resulted negligible 
and is not included in the plots. 

For the bottom two of the four plots panels, the curves 
are divided by the actual value of the parameter in order 
to express the relative uncertainty. The general trend 
is that the CRLB fails to describe correctly the error 
for SNRs lower than 20. For the t and the the M this 
regime is of physical interest while for the symmetrized 
mass ratio and the phase the CRLB already predicts very 
large uncertainties. It is also worth noticing that the SNR 
at which the second order covariance matrix contribution 
becomes comparable to the CRLB is much larger when 
the full parameter space is involved. For t for example 
in the NS-NS case the two are the same a.t p — 2 while 
for the full parameter space they equate around p = 15. 
This results appear to indicate that also the timing based 
directional reconstruction accuracy, is worse when the 
other physical parameters of the waveform are unknown. 
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Figure 3. NS-NS signal in initial LIGO noise . The dotted 
line is the CRLB. The dashed line is the square root of the 
second order covariance matrix, and the continuous line is the 
square root of the sum of the diagonal elements of the FIM 
and of the second order covariance matrix. In the last two 
panels the errors are divided by the value of the parameter. 
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Figure 4. NS-BH signal in initial LIGO noise. The dotted 
line is the CRLB. The dashed line is the square root of the 
second order covariance matrix, and the continuous line is the 
square root of the sum of the diagonal elements of the FIM 
and of the second order covariance matrix. In the last two 
panels the errors are divided by the value of the parameter. 



Figure 5. BH-BH signal in initial LIGO noise. The dotted line 
is the CRLB. The dashed line is the square root of the second 
order covariance matrix, and the continuous line is the square 
root of the sum of the diagonal elements of the FIM and of 
the second order covariance matrix. In the last two panels 
the errors are divided by the value of the parameter. 
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Figure 6. NS-NS signal with Advanced LIGO noise . The 

dotted line is the CRLB. The dashed line is the square root 
of the second order covariance matrix, and the continuous 
line is the square root of the sum of the diagonal elements 
of the FIM and of the second order covariance matrix. In 
the last two panels the errors are divided by the value of the 
parameter. 




Figure 7. NS-BH signal with advanced LIGO noise. The 

dotted line is the CRLB. The dashed line is the square root 
of the second order covariance matrix, and the continuous 
line is the square root of the sum of the diagonal elements 
of the FIM and of the second order covariance matrix. In 
the last two panels the errors are divided by the value of the 
parameter. 
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Figure 8. BH-BH signal with advanced LIGO noise. The 

dotted line is the CRLB. The dashed hne is the square root 
of the second order covariance matrix, and the continuous 
line is the square root of the sum of the diagonal elements 
of the FIM and of the second order covariance matrix. In 
the last two panels the errors are divided by the value of the 
parameter. 
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Figure 9. The simulations of T. Cokelaer (first and thrid 
panel from top), compared with our result (see the text for a 
discussion) . 
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NS-NS SNR 20 10 6 



Cokelaer 


25% 


200% 


700% 


Z.V.M. 


8.6% 


131% 


174% 


BH-BH SNR 


20 


10 


6 


Cokelaer 


10% 


150% 


230% 


Z.V.M. 


2.9% 


111% 


129% 



Table I. The discrepancies between the CRLB error predic- 
tion and Montecarlo simulations are presented above the dis- 
crepancies observed using the first and second order covari- 
ance matrix. The discrepancies are presented for three values 
of the SNR: 20,10 and 6. 



It interesting to compare our results with Monte Carlo 
shuulations for 3.5PN waveforms, (^|59|). 

In fig. [9]we reproduce the plots of ([S^) obtained for a 
NS-NS system (first panel from top) and for a low mass 
BH-BH system - lOM© - (third panel from to) when the 
initial LIGO noise is used. 

The results we obtain for the same systems are shown 
in the second and fourth panels from the top of fig. [9j 
where the CRLB and the square root of the sum of the 
inverse FIM and the second order covariance are plot- 
ted. The plots show that the use of the second order 
covariance predicts correctly the SNR regime where the 
CRLB fails to describe the MLE accuracy. The explicit 
percentage discrepancies presented for SNR equal to 20, 
10 and 6 in table |T] seem to indicate that higher orders of 
the expansions might be necessary to fully reproduce the 
error derived with the Monte Carlo simulations. 



VII. CONCLUSIONS 

In this paper we applied a recently derived statisti- 
cal methodology to gravitational waves generated by the 
inspiral phase of binary mergers and for noise spectral 
densities of gravitational wave interferometers. Explic- 
itly we computed the first two orders of MLE expansions 
of bias and covariance matrix to evaluate MLE uncer- 
tainties. We also compared the improved error estimate 
with existing numerical estimates. The value of the sec- 
ond order of the variance expansions allows to get error 
predictions closer to what is observed in numerical simu- 
lations than the inverse of the FIM. The condition where 
the second order covariance is negligible with respect to 
the first order predicts correctly the necessary SNR to 
approximate the error with the CRLB and provides new 
insight on the relationship between waveform properties 
SNRs and estimation errors. Future applications include 
IMR waveforms, network of detectors and source location 
estimation. 
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IX. APPENDIX A: DERIVATION OF THE 
EXPANSIONS 

Analytic expressions for the moments of a (MLE) are 
often difficult to obtain given a non-linear data model. 
However, it is known from ([1? ) that likelihood expan- 
sions can be used to obtain approximate expressions for 
the moments of a MLE in terms of series expansions in 
inverse powers of the sample size (n). These expansions 
are valid in small samples or SNR where the MLE may 
have significant bias and may not attain minimum vari- 
ance. 

An expressions for the second covariance equivalent to 
the one used here, was first given in a longer form (about 
a factor 2), and with the prescription of isolating the dif- 
ferent asymptotic orders inside the tensors, in ( JQJ). The 
expression presented here were first derived in an unpub- 
lished MIT technical report ([35]) by two of the authors 
of this paper, where (a) the Bartlett identities (|,12, ,^13j) 
were used to simplify the expression of the second or- 
der variance, and derive the second order bias, (b) the 
prescription on the tensors no longer needed to be imple- 
mented. The final form of the second order covariance 
has already been published without proof, in two papers 
([37], [35]), were the first and third author of this paper 
are coauthors, involving MLE of source and environment 
parameters that use acoustic propagation within shallow 
water ocean wave guide. 

In this section, we derive the second order bias for 
a multivariate MLE and we introduce a chain rule 
that allows the derivation of the of the second order 
covariance matrix from the second order bias. The 
explanation follows closely [39]. The derivation of 
the bias is performed in two steps: first we derive the 
expansions in the non physical scenario of n statistically 
independent identical measurements and then set n=l 
for the case of interest of this paper. 

The derivation follows the approach of ([I^) for likeli- 
hood expansions, (|21j) for the large sample approxima- 
tions and ([3], [TT]) for the notation of the observed like- 
lihood expansions. For the small sample case of interest, 
asymptotic orders generated from the likelihood expan- 
sion may contain different powers of ^ , and the contribu- 
tions to a single power of ^ may have to be collected from 
different asymptotic orders. ([26]) extending the work of 
([17j). confronted the equivalent problem of separating 
the powers in n within the expectations of products of 
linear forms and arrived at expressions for the second or- 
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der bias of the MLE by embedding the derivation in a 
discrete probabiHty scheme. Some apphcations of their 
results are given in ([15]), ([55]), ([Tn]),([Il]),([in]). We 
present here an independent derivation for the second or- 
der bias that is vahd for general discrete or continuous 
random variables. 

Let us consider a set of n independent and identically 
distributed experimental data vectors = {xn^ ..,XiN}, 
where N is the dimension of every vector. We assume 
that these data are described by a probability density 
Pix:"^) = nip(xj,i9) that depends on a D-dimcnsional 
parameter vector i? = According to ([H]), 



where (z9 — i?)'" for r = 1, .., Z? are the components of the 
observed error. We will use the notation 

ia2..as •• 

where Ha^a2..a, = ^ama-.a, " Waiaa-.Q. and -Vab is the 



we suppose that the MLE = {t9i, .., i^^)} of -d is 
given by a stationary point of the likelihood function 
K-a^iQ — l'n{P{x,'d)) = Xir=i ^"(pfei'^)) ^'^^^ respect 
to the components of "ff 

'.fel) = ^^l,.g = r^l,..,D. (9.1) 

The first step in deriving the likelihood expansion, if 
Ir (x, 1?) can be expanded as a Taylor series in the com- 
ponents of the observed error — i?, consists of writing 
lr{x,d) as 



(9.2) 

I 

information matrix iab = —Vab — ~E[lab\ — E[lalb]- In- 
troducing as the inverse of the matrix whose elements 
are given by jab = —lab, equation ( |9.2[ ) can be rearranged 
to solve for {dr - ^r) = (d-^Y 



= lr{x,d)=^lrix, {i-l)+l)=lrix,l)+lrsix,l)id~^y + -lrtu{x,l)i'^-^Y{-^-^T + --- 



{d-dy = rh + ^rht^{d- ^)\d- ^y + -fHstuvid-^y{d-^y{d- i)y 

+ :^ristuvU^-i)y{d~dy{d~dy{d-dr + ... . (9.3) 



Finally we can iterate equation (9.3) with respect to 
the components of the observed error and expand j"^ 
in terms of the information matrix inverse — {i^^)ab 
(0, page 149) 



(9.4) 



From (9.3) and (9.4) the terms that contribute to each 



asymptotic order of the observed error can be obtained. 
However, in order to isolate the asymptotic orders nec- 
essary to compute the second order bias we ha ve ch osen 
only a finite number of terms within equations (9.3) and 



(9.4). This choice can be made by (a) observing that 



is proportional to n 

mi 1 d'ln{p{xi,l)) 

where the value of i is irrelevant because all the vector 
data have the same distribution, and (&) using the large 
sample approximation expressed by ([21], page 221) 



In 



E\L 



n2 



(9.6) 



and proved for discrete and continuous random variables 
in the next paragraph. 

Equation (9.6) indicates that the expectation of a 



product of c generic factors Ha^a^-.a^ is a polynomial of 
integer powers of n, where the highest power of n is the 
largest integer less than or equal to | (we use the nota- 
tion int{^)). For example, U{Qia2..a J{6ib2..hp}{cic2..c,} is 
proportional to n. 



The proof of the large sample approximation (9.6) is 



an extension to continuous random variables of the anal- 
ysis performed in ((27j,page 36), for discrete random vari- 
ables. 

To prove equation (|9.6| we show that the quantities 



(9.7) 



obtained as the expectation of products of m functions 
Ha. = X)"=i l^'i.a^ are polynomials in integer powers of n 
of order less than or equal to ^. 



Here, /i, 



I, 



< li^a- >. The subscripts that 



''2, a 

appear in Ha represent collections of indexes as intro- 
duced in equation (9.6). 

The factors Ha^Ha^---Ha appearing in the expectation 
can be seen as the first m terms of a succession Ha . , 
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where we choose arbitrary Ha for i > m. Using such a where we assume that all the expectations exist, 
succession we can define the quantity 

oo CO n n Moving all the sums outside of the expectations and 

F{y) — ''^^y'^ Ps{n) = ^^J/^'-^K^^ ^js,aj] redefining the notation in the sums, we can obtain 

s=0 s=0 ji = l js=l 



Fiy) 



oo 
s=0 



,i+..+ai. 



, +..+a„ „ =s 



(9.8) 



where the quantity ap^a^ can be equal to or to 1 to 
indicate that the quantity hp^a is present or not in the 
product inside the expectation. The summation over all 
ai,a. is intended to account for all possible configurations 
of the indexes (choosing some equal to zero and the oth- 
ers equal to one) with the constraint that the number of 
factors within the expectation is equal to s. We can now 
group together terms (/ii,Q .)°' °^ for different i and de- 



fine o-i,a- = Sii where ~ ^^is manner 

we can also factorize the expectations in (9.8) as expec- 



tations for different data vectors. By making use of the 
statistical independence of the different data vectors, and 
defining A(ai,Q. ,...,ai,Q. ) = E[hi^a. ■■■ht,a 

— 31 — 3p — Ji — : 

can be rewritten as follows: 



(9.8) 



S— Si—0 a-i^ai + ■■ -+131,0^ —Si S2—0 a2^ai + .--4-a2,a^ —S2 

• • • ^ ^ ^ ^ ^(^^71, 7 ■ ■ • ; ^n,a_^ ) 

Sn — S — Si—S2 — --. — Sn-i a^^^i + + S„ 



"(E^'' E A(ai,„^, . . . ,ai,aj 

=0 E^=i ai,Qj.=si 



(9.9) 



We observe that, when the expectations contain only 
one factor, for si 1, we have 

E[hiJ^^-^ihi^^^^-^2...] = E[hi^a^] = A(ai,a^.) = for any 
j. As a consequence, the Taylor expansion of f{y) in y 
can be written as f{y) = cq + C2y^ + c^y^ + ... 
Eventually, it is the absence of Ciy in this expansion that 
allows us to explain the properties of the polynomials in 
equation (9.7) and finally the large sample approxima- 
tion (9.6 1. To accomplish this last step we explain how 



the polynomials P,„(rt) are related to the coefficients 
Let us consider the contribution to Pm{n) that originates 
from the product of k factors Ci-^y^^ , .., Ci^y''' with the 
constraints ii > 1 and ii + .. + ik = m. The number of 
ways these products can be formed for an assigned set of 
ii, ..,ik is proportional to (^) ~ n*^ -I- lower powers in n. 
Moreover, the contributions to Pm{n) are formed by 
an integer number of factors less than or equal to 
Y because ci — 0. This property limits the highest 
power in n contained in P„i(n) with the largest integer 



number smaller than or equal to ^. Equation (9.7) 
is then proved as is the large sample approximation (|9.(j[). 



Each asymptotic order m of the r*'* component of the 
observed error is denoted by (-i? — i?)''[to] where the index 



is given by m = ^i— and s is a natural number. For 



— ^+1 

2 

example, the asymptotic order for m = 1 is given by 
1 



(^-^)'-[l] = ^v-'^Hjt 



(9.10) 



where we have adopted the lifted indexes notation 
ya...z _ _^i^Sy^ xhe asymptotic orders of the bias 
are then given as the expectation of the asymptotic or- 
ders of the observed error 



(9.11) 



The asymptotic order b['d^,m\ contains different powers 
of ^ as we discuss in this paragraph. 



It follows from equations (9.3) to (9.6) that b['d'^,m] is 
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the sum of terms having the structure 




where a, 6, and c are 



the numbers of factors in the three products satisfying 
a — 6 — I = m. Different powers of ^ can be gener- 
ated because E[H(^ 'f .. can contain ah the integer 

c 

powers of n equal to or less than However 



can 



from the fact that no power higher than n 
be generated follows that an asymptotic order m will 
never generate powers in the sample size with p 
smaller than m. It still remains to prove which is 
the largest power in ^ contained in the asymptotic 
order m. We show that in the following. Since the 
largest range of powers of ^ is allowed when & = 0, we 
study these terms to evaluate the range of powers of ^ 
contained in b['d^,m\. The structure of equation (9.3) 



implies that its iteration with respect to the observed 
error components generates an equal number of j""^ 
and -ff ( ) (we recall that ^( ) = i7( ) + ^^(.)). Similarly, 
the number of i^'^ generated from the expansion of j"'' 



given in equation (9.4) is equal to the number of Hab 
plus 1. These two observations imply that the terms 
where 6 = also verify a = c = 2m. As a consequence, 
the powers of ^ that 6['!9'",m] contains can range from 
the smallest integer number larger than or equal to m, 
up to 2m — 1. 

The analysis above implies that in order to compute 
the contributions to ^ and to ^ in the bias expansion 

that we denote with &[i?'',l] and 6[i?'",2], it is necessary 
to obtain the first three asymptotic orders of the bias, 
b[-d^, 1], b[d^, |], and 5[??'',2]. In the explicit expressions 
below, to condense the notation, we introduce the 
quantities = i"'^..^'^ , so that 



I 2 '-'ap-y ^st,u,v ' r2 ^stv,u.w > ^st^vu^w ~i~ ^-^ ^st,w 

I jra.st.bc.de^ , , ^ jrj.kl.rap.qz ^, , r ^ ^,roLB'y5 , 

+ -'''■ W{af,}{cd}{et}{s} + g-' "{j;p4{fe}{™}{-3} + + 



(9.12) 
(9.13) 



(9.14) 



The first order of the bias 1] can than be obtained 
by substituting the explicit expressions for the tensors 
in 6[i?'', 1]. The second order 2] takes contributions 
from bl^"^ , |] and 6[i9'',2]. However, while bld"^ , |] gener- 
ates only n^^ contributions, b['d'^ ,2] generates n^"^ and 
n~^ contributions. Consequently, to collect all and only 
the contributions to n~^, it is necessary to extract the 
n^^ component from 6[i?'",2] and add it to 6[i9'', |]. The 
extraction can be done by introducing into 6[i?'',2] only 
the highest powers of n of the tensors. 

The first two orders of the bias for the MLE of the r 



component of the parameter vector i9 then become (2.2) 



and (2.3) 



Starting from (2.3) the form for the second order 



covariance matrix can be obtained using the following 
theorem: 

Theorem: Let g{.) he an invertible and differentiable 
function and let ^ = g{^) be a set of parameters depen- 
dent on ^. Let 6(??r) o.nd 6(^r) be the biases relative 
to the estimation of the components of ^ and ^ which 
can be computed as power series of ^, as explained in 
Section 2. The terms of the expansion in powers of ^ 
forb{^r) can be obtained from the terms of the expansion 
for 6('i?,.) by replacing the derivatives with respect to 
the components of with derivatives with respect to 
the components of ^. The explicit dependence on ^ can 
be removed at the end of a derivation by means of the 
Jacobian matrix J defined by 



d 



a-d. 



d 
d-ds 



J 

■Jms 
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To prove the chain-rule theorem given in section [IT], we 
analyze how the derivation of the expansion in powers of 
^ for b{'dr) r = 1, .., D, which is described there, changes 
if we are interested in the expansion for f — ^, ■■, D 
where ^ is a vector of parameters in an invertible relation- 
ship with t- fit); t= f'^ii)- The starting point 
for the derivation of the expansion of fe(^r) t = ^, D is 
the stationarity condition 







1, 



,£1.(9.15) 



that can be obtained from equation (9.1) by replac- 
ing only the derivatives gg- forr = 1,..,D with 

g|- for r = 1, .., D, because 'd = t implies ^ — f(d) ~ ^. 
The subsequent steps can then be obtained from equa- 



tions (9.2) up to (2.3) by replacing the derivatives in the 
same way. Since = ^afl ^ "^"'^afl' where the 
components of the Jacobian matrix Jms behave like con- 
stants in the expectations, the substitution of the deriva- 
tives can also be done only in the final result (for example, 
in the orders given in equations (9.14) and (2.3)). The 
expectations contained in the expansion of 6(i?7) in pow- 
ers of ^ can also be computed before the substitution of 
the derivatives if the likelihood function dependence on 
the parameters is expressed in terms of auxiliary func- 
tions. Examples of auxiliary functions are griil) = i^r 
for a general parametric dependence and gi{i)) = /i(i9) 
52(1^) = ('■(i^) for a scalar Gaussian distribution. By 
means of these auxiliary functions, the derivatives 

and g|- become g|- and g|- . As a consequence 
the orders of the expansion for 6(^r) can be found from 
the orders for the expansion of 5(t?r) implementing, in 
the result of the expectations, the substitutions 



8''+''+-+'- grail) 9'^+*^ + -+''°.9m(/-^(e)) 



{d'di)^K.{d'&D)'" 



(aci)^i..(a^z?^- 



(9.16) 



The converse of the chain-rule, in which higher moments 
of the observed error are used to compute lower moments, 
is not possible. We can observe, for example, that the ex- 
pansion of a general moment of an MLE does not always 
contain all the powers of ^. The lowest order present 
in the expansion of the m order moment is given by the 
largest integer number smaller than or equal to ^ (we 
use the notation int{^^)): 



it(f ) 



higher powers of — . 

n 



The consequence of this observation is that only the bias 
and the error correlation matrix may contain first order 



terms. For this reason an inverse chain-rule would have 
to generate non-zero orders of lower moments expansions 
from the corresponding orders of the higher moments ex- 
pansions that are zero for powers of - lower than . J", m-, . 

Let us consider how the chain-rule makes it possible 
to compute the expansion of the error correlation ma- 
trix ■~{d) defined by = E[0_ - d)Q_ - -d)^] and the 
covariance matrix. Using the invariance property of the 
MLE, 



S(i?) = b[{t£) ~ d{b{l))^ - b(l)if 



(9.17) 



where (.)"^ is the transpose of a vector and b{&d^) is a 
matrix whose components are the bias of the product of 
two components of ^. Once S(;^) and the bias are known, 
then the covariance matrix C('!9) can also be computed 
by means of 



c{d) = E{d) - b{d)b{dy' 



(9.18) 



To compute the right hand side of equation ( 9.17[ ) we 
express it in terms of the components, obtaining 

E,,id) = bi^j) - d.bidj) - b{d,)d,. (9.19) 

It is important to realize that knowledge of 6(i?r) is suf- 
ficient because the expansion of 6(i?ii?j) can be derived 
from it by means of the theorem given above. In fact, if 
we choose 



(9.20) 



as a new set of parameters, b{'di'dj) becomes b{£,i). How- 
ever, is necessary to insure that the relationship between 
z9 and ^ is invertible. This condition holds if both 7^ 
where i = 1, .., Z? and the sign of the components of d is 
known, for example, by means of constraints on the data 
model. 



In a scalar estimation scenario, equation (9.17) becomes 
simply 

E{d) = E[{d~ ^ b{^) - 2db{d) (9.21) 

and the variance C{d) = 2(1?) — b{d)'^. In this case b{{}'^) 
can be derived from b{d) if we choose f = i?^ as the new 
parameter and the Jacobian matrix becomes ^ because 
d_ _ M jL — J_JL 

a? ~ ~ 2i? d€ ■ 

A useful simplification of the algebra of the chain-rule 
in the derivation of second order moments is described 
in the following two paragraphs. The chain-rule and the 
subsequent conversion of the derivatives require the sub- 
stitutions 



gr, 



r I o 



M) 



J2 



(9.22) 
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From the right hand side of ( |9.22 1 it is clear that the 
derivatives will also be applied to the Jacobian matrix, 
thereby generating m! terms for every derivative of 
order m. However, it can be observed that the terms 
contributing to the final result are only a small subset. In 
fact, among all the terms generated from the conversion 
of the derivatives in the bias expansion, only those where 
a first derivative of the Jacobian matrix appears must 
be considered. For example, we can consider the term 

u 



ra,sb,vc 



jI^'^''^'^'^-^VdefVs,t,u that comes from 
equation (2.3), which must be used to derive the second 
order In this case, we need to consider only the 

3 terms in which one of the derivatives represented by 
a,b,c operates on a Jacobian matrix and d,e,f operate 
on the likelihood function, plus the 3 terms where the 
role of a, b, c and d, e, / are inverted. In general, the 
total number of relevant terms among those generated 
in every derivative is equal to or less than the order of 
the derivative m. 



The detailed analysis of equation (9.19) reveals that 
the terms generated in hiasl^i-dj) can be divided into 
three groups: (a) the terms where no derivative of the 
Jacobian matrix appears are equal to {{J~^)'^)isbias{'ds) 
(we show in Example 1 that {{J~^)'^)isbias{'ds) cancels 
with -~'dib('dj) — b{'di)i!)j after its introduction in equation 
( |9.19 )); (b) the terms where only one derivative of the 
Jacobian matrix appears give the error correlation 
matrix; and (c) the terms that contain higher derivatives 
of the Jacobian matrix or more than one first derivative 
of the Jacobian matrix in the same term, summed 
together, give zero. 

To clarify the use of the chain-rule and the alge- 
braic issues discussed above, we present two examples. 
In Example 1 we use the first order term of the bias 
in a general vector estimate to derive the first order 
term of It is useful to recall that the expansion 

of S(i9) and the expansion of C{d) can be different 
only from the second order on, so this example also 
describes a^derivation of the first order of the covariance 
matrix C(i?). Following the same approach, the second 
order term of the error correlation matrix expansion can 



be derived from equation (2.3 1 and the second order 



covariance matrix can also be derived if we also use (2.2) 



and (9.181 



In Example 2 we illustrate the way the chain-rule can 
still be used, if the available expression of the bias 
expansion is explicitly computed for a particular data 
model. In particular, we derive the second order mean 
square error and the second order variance from the 
second order bias in two scalar Gaussian models. In Ex- 
ample 2 we also illustrate the simplification introduced 
above for the algebra involved in the conversion of the 
derivatives. 



Example 1 

Using the Bartlett identities ([ 



]), we rewrite the 



first order bias, given by equation (9.12), as 



(9.23) 



From equation (9.231, b{-di'dj) = b{^i) can be computed 
by means of the chain-rule by replacing the derivatives 
with respect to the components of z9 with derivatives with 
respect to the components of ^, where ^ is given in equa- 
tion (9.20) and using the corresponding Jacobian matrix. 



The tensors in equation (9.231 become 



sp 



JjaJll3JpfVa,j3,j{'&) 



sp 
-1\ 



where we have specified in the bracket beside the tensors 
the dependence on the parameter sets. Inserting these 
expressions in equations (9.23) and observing that equa- 



tion (9.19) can be expressed in the form 

S(|),, = b{^/) - {{J-Yhsb{ds) (9.24) 

the first order term of the error correlation matrix can 
then be obtained as 



2" 5i?/3 ~ 2' 
where we have introduced the tensor 



■a/3 



7a0 



(9.25) 



7a/3 



1 iia = if5 = joTa~j/3~i 
otherwise. 



Example 2 

In this example we determine the second order variance 
and mean square error from the second order bias for two 
cases of parameter dependence for the scalar Gaussian 
density 



p(x,i9) 



1 



(27rc)"/2 



exp 



1 " 



(xi - fi)'^ 



(9.26) 



as a direct application of equation (9.21 1 



In the case where c does not depend on the param- 
eters (^ = 0) the second order bias can be derived 
using the scalar version of equation (2.3) and of the 
tensors v defined in equation (9.5). For this parameter 



dependence, the asymptotic order for to = | (equation 



(9.13)) is zero, and the second order bias can be di- 
rectly obtained also from equation ( 9.14[ ) The result is 



b[M = (f^ir/^XA)' - I'M (A)' - fKi^n Applying 
the chain-rule, the second order mean square error for 
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the estimation of § becomes S(i?)[^] = — j^, 

where the full conversions of the derivatives are given by 
f, . _dlL - JL ■ r, ^ - 1 

ci3,, 1 , ir. -i,; > _ 

(ai52)4 — 
By means of these con- 



"fi' _ 3 'A' I 15/i 



_3/i_ , 



versions and the expression of the second order bias, it 
can be observed that among the 18 terms that are in 
principle generated by the chain-rule, only 6 contribute 
to the second order mean square error. 

In the following we show that the expansion in the 
inverse of the sample size is equivalent to an expansion 
in 1/SNR. The derivation of the asymptotic orders in 
- would be the same for an expansion in any quantity 



if (9.5) and (9.6) can be derived for a certain quantity 7 



instead of the sample size n. In this section we illustrate 
indeed that this is the case where the signal to noise 
ratio for a set of scalar data distributed according to 
a deterministic Gaussian PDF takes the role of the 
sample size. The probability distribution and parameter 
dependent part of the likelihood function are given by: 



1 



(9.27) 



(9.28) 



We also define the signal to noise ratio for this exam- 
ple as 7 = /x^/tT^ following standard practice for scalar 
deterministic Gaussian data. We can obtain: 



iJ„(^) = 2^„(^)Si(^i^iM (9.29) 



(l?) = > C{p) + W 5 



(9.30) 



p r- r~ 



where a or p denote the order of the derivative or an 
arbitrary set of derivatives. 

^ ^ n 

From (|9.31| it becomes obvious that Va is proportional E[Ha, ■ ■ ■ H^J = 2^^^^ ... Map ^^[JZ^^* ~ 



i=i 



to not only the sample size, n, but also to the 7. The ^ 

term inside the square brackets is simply a sum over nor- _ 2P(^^f^l) ^ ^°'p ) (^2)' ^ 

malized derivatives of the mean and contains information ^ a^P 25 (|)! 

about the shape of the signal. In the above equations we 2 

use: =(^)2...( )...(^)n2 



^ E[Ha,...Ha^]^{n-iY^ (9.33) 

It is therefore found that E[Hat---Hap] is propor- 
(o„\\ tional to (7)**^^ and we have proved the analogy between 

E[{x - m(i9))^"] - (9-32) sample size and 7. 

Note that for deterministic Gaussian data, the non- 
integer asymptotic orders of the bias are zero and the 
integer orders are equal to 61, 62 etc. This is sufficient to 
prove that the orders in 1/n of the bias expansion are also 
orders in I/7. This also holds for the variance expansion. 
The next step is to determine the dependence of A similar, although longer proof can be written for the 
E[Hat ■ ■ ■ Ha ] on 7. This is shown below: SNR definition provided in|4.5[ 
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